NicheNet on MM cells (NK cells) using ligands found by LIANA:

library(Seurat)
## Attaching SeuratObject
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.1.8
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(magrittr)
## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(liana)
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(BiocParallel)
library(nichenetr)

Load Rds with NK cells, liana results and filter the interactions by ICR

dat <- readRDS("Data/MM_atlas.Rds")# scRNA-Seq atlas
receiver = c("NK exhausted", "NK resident")
## labels
library(readr)
all_labels_df <- read_csv("all_labels_df.csv", 
    col_types = cols(...1 = col_skip()), 
    trim_ws = FALSE)
## New names:
## • `` -> `...1`
all_labels_df <- as.data.frame(all_labels_df)
rownames(all_labels_df) <- all_labels_df[,"keyName"]
all_labels_df$keyName <- NULL
colnames(all_labels_df) <- "labels_interactions"
all_labels_df$labels_interactions <- as.factor(all_labels_df$labels_interactions)

dat <- AddMetaData(dat, all_labels_df)

dat <- SetIdent(dat, value = "labels_interactions")
ICR = c("ADCY7|ADGRE5|ADGRG1|BTN3A2|CD244|CD47|CD53|CD81|CD96|CLEC2D|CXCR4|DIP2A|FAS|HAVCR2|HLA-DPB1|IGF2R|KIR2DL1|KIR2DL3|KLRB1|KLRC1|KLRC1_KLRD1|KLRF1|KLRG1|LAG3|LAIR1|LILRB1|PTPRC|RACK1|SIGIRR|TIGIT|TNFRSF14")

interactions <- read_csv("all_interactions_defaults.csv", col_types = cols(...1 = col_skip())) %>%
  filter(target %in% receiver) %>%
  filter(grepl(ICR,receptor.complex)) %>%
  filter(!grepl("KLRC2|KLRC3", receptor.complex))
## New names:
## • `` -> `...1`

Load models for NicheNet:

ligand_target_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_target_matrix_nsga2r_final.rds"))

ligand_tf_matrix = readRDS(url("https://zenodo.org/record/7074291/files/ligand_tf_matrix_nsga2r_final.rds"))

lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))

sig_network = readRDS(url("https://zenodo.org/record/7074291/files/signaling_network_human_21122021.rds"))

gr_network = readRDS(url("https://zenodo.org/record/7074291/files/gr_network_human_21122021.rds"))

weighted_networks = readRDS(url("https://zenodo.org/record/7074291/files/weighted_networks_nsga2r_final.rds"))

weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))
# Select ligands actively targetting ICR
ligands = unique(sort(interactions$ligand.complex))
ligands = ligands[ligands %in% colnames(ligand_target_matrix)]
ligands
##  [1] "B2M"      "CALM1"    "CCN2"     "CD22"     "CD40LG"   "CD48"    
##  [7] "CD55"     "CD99"     "CDH1"     "CEACAM1"  "CLEC2B"   "CLEC2D"  
## [13] "CXCL12"   "FAM3C"    "GNAS"     "GZMB"     "HLA-A"    "HLA-B"   
## [19] "HLA-C"    "HLA-DPB1" "HLA-DQA1" "HLA-DQA2" "HLA-DQB1" "HLA-DRA" 
## [25] "HLA-DRB1" "HLA-DRB5" "HLA-E"    "HMGB1"    "IL1B"     "LGALS1"  
## [31] "LGALS3"   "LGALS9"   "LILRB4"   "LTF"      "MIF"      "NECTIN1" 
## [37] "NECTIN2"  "NRG1"     "SIRPA"    "TGFB1"    "TGM2"     "THBS1"   
## [43] "TNFSF13"  "TNFSF13B"

Define baseline gene expression for NK cells:

dat_nk = subset(dat, subset = labels_interactions %in% receiver)

I will consider a gene expressed if it is present in at least 10% of cells for each celltype as suggested in NicheNet Documentation

## receiver
expressed_genes_receiver = get_expressed_genes(receiver, dat_nk, pct = 0.10)

background_genes = expressed_genes_receiver %>% .[. %in% rownames(ligand_target_matrix)]

Select gene set of interest

As target genes I am going to use the genes included in those regulons that contains ICRs mediators from the TF analysis. This is aimed to integrate the two endpoints, signalling and regulatory network, to buil the signalling cascade of ICRs:

geneset_oi <- read_tsv("Data/all_ICR-mediators_target_genes.csv", col_types = cols(), col_names = "gene") %>%
  pull(gene) %>%
    .[. %in% rownames(ligand_target_matrix)]

Run NicheNet

nichenet_activities <- predict_ligand_activities(
  geneset = geneset_oi,
  background_expressed_genes = background_genes,
  ligand_target_matrix = ligand_target_matrix, potential_ligands = ligands
)

Visualize ranked interactions by LIANA and NicheNet:

# prepare data for visualization
vis_liana_nichenet <- interactions %>%
  inner_join(nichenet_activities, by = c("ligand.complex" = "test_ligand")) %>%
  arrange(pearson) %>%
  mutate(ligand = fct_inorder(ligand.complex))

# prepare NicheNet figure
nichenet_scores_plot <- vis_liana_nichenet %>%
  group_by(ligand) %>%
  summarize(pearson = mean(pearson)) %>%
  ggplot(aes(y = ligand, x = pearson)) +
  geom_bar(stat = "identity") +
  ggtitle("NicheNet") +
  xlab("Pearson's score") +
  theme_cowplot() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.line.y = element_line(color = "white"),
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))

# prepare LIANA figure
liana_receptor_heatmap <- vis_liana_nichenet %>%
  ggplot(aes(y = ligand, x = receptor.complex, fill = aggregate_rank)) +
  geom_tile() +
  theme_cowplot() +
  ggtitle("LIANA") +
  ylab("Ligand") + xlab("Receptor") +
  theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1),
        plot.title = element_text(hjust = 0.5),
        panel.grid.major = element_line(colour = "gray", linetype = 2),
        legend.position = "left")

# combine plots
plot_grid(liana_receptor_heatmap, nichenet_scores_plot,
          align = "h", nrow = 1, rel_widths = c(0.8,0.3))

Rank ligands by pearsons score

nichenet_activities %>% arrange(-pearson) 
## # A tibble: 44 × 5
##    test_ligand auroc   aupr aupr_corrected pearson
##    <chr>       <dbl>  <dbl>          <dbl>   <dbl>
##  1 THBS1       0.657 0.0634        0.0566   0.134 
##  2 LGALS3      0.589 0.0546        0.0478   0.0976
##  3 CD40LG      0.707 0.0309        0.0241   0.0876
##  4 CCN2        0.630 0.0530        0.0462   0.0717
##  5 TNFSF13     0.629 0.0201        0.0133   0.0708
##  6 HLA-DRB5    0.634 0.0260        0.0192   0.0666
##  7 CD22        0.666 0.0191        0.0123   0.0629
##  8 HLA-C       0.681 0.0175        0.0107   0.0601
##  9 LGALS9      0.666 0.0164        0.00958  0.0552
## 10 SIRPA       0.667 0.0260        0.0192   0.0544
## # … with 34 more rows
best_upstream_ligands = nichenet_activities %>% arrange(-pearson) %>% pull(test_ligand)

Predict target genes for each ligand

active_ligand_target_links_df = best_upstream_ligands %>% lapply(get_weighted_ligand_target_links,geneset = geneset_oi, ligand_target_matrix = ligand_target_matrix, n = 250) %>% bind_rows()

active_ligand_target_links_df = na.omit(active_ligand_target_links_df)

active_ligand_target_links = prepare_ligand_target_visualization(ligand_target_df = active_ligand_target_links_df, ligand_target_matrix = ligand_target_matrix, cutoff = 0.25)

nrow(active_ligand_target_links_df)
## [1] 100
head(active_ligand_target_links_df)
## # A tibble: 6 × 3
##   ligand target weight
##   <chr>  <chr>   <dbl>
## 1 THBS1  BCL2L1 0.0152
## 2 THBS1  PECAM1 0.0668
## 3 THBS1  SOCS3  0.0117
## 4 LGALS3 BCL2L1 0.0102
## 5 LGALS3 SOCS3  0.0649
## 6 CD40LG BCL2L1 0.0252
order_ligands = intersect(best_upstream_ligands, colnames(active_ligand_target_links)) %>% rev()
order_targets = active_ligand_target_links_df$target %>% unique()
vis_ligand_target = active_ligand_target_links[order_targets,order_ligands] %>% t()

p_ligand_target_network = vis_ligand_target %>% make_heatmap_ggplot("Microenvironment ligands","Inhibitory target genes", color = "purple",legend_position = "top", x_axis_position = "top",legend_title = "Regulatory potential") + scale_fill_gradient2(low = "whitesmoke",  high = "purple", breaks = c(0,0.005,0.01)) + theme(axis.text.x = element_text(face = "italic"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
p_ligand_target_network

Infering ligand-receptor activities

# get the ligand-receptor network of the top-ranked ligands
receptors = lr_network %>% pull(to) %>% unique()
expressed_receptors = intersect(receptors,background_genes)
lr_network_top = lr_network %>% filter(from %in% best_upstream_ligands & to %in% expressed_receptors) %>% distinct(from,to)
best_upstream_receptors = lr_network_top %>% pull(to) %>% unique()

# get the weights of the ligand-receptor interactions as used in the NicheNet model
lr_network_top_df = weighted_networks$lr_sig %>% filter(from %in% best_upstream_ligands & to %in% best_upstream_receptors)

# convert to a matrix
lr_network_top_df = lr_network_top_df %>% spread("from","weight",fill = 0)
lr_network_top_matrix = lr_network_top_df %>% select(-to) %>% as.matrix() %>% magrittr::set_rownames(lr_network_top_df$to)

# perform hierarchical clustering to order the ligands and receptors
dist_receptors = dist(lr_network_top_matrix, method = "binary")
hclust_receptors = hclust(dist_receptors, method = "ward.D2")
order_receptors = hclust_receptors$labels[hclust_receptors$order]

dist_ligands = dist(lr_network_top_matrix %>% t(), method = "binary")
hclust_ligands = hclust(dist_ligands, method = "ward.D2")
order_ligands_receptor = hclust_ligands$labels[hclust_ligands$order]
vis_ligand_receptor_network = lr_network_top_matrix[order_receptors, order_ligands_receptor]
p_ligand_receptor_network = vis_ligand_receptor_network %>% t() %>% make_heatmap_ggplot("Microenvironment ligands","Inhibitory receptors", color = "mediumvioletred", x_axis_position = "top",legend_title = "Prior interaction potential")
p_ligand_receptor_network

Network

Complete network

active_signaling_network = get_ligand_signaling_path(ligand_tf_matrix = ligand_tf_matrix, ligands_all = ligands, targets_all = geneset_oi, weighted_networks = weighted_networks)
data_source_network = infer_supporting_datasources(signaling_graph_list = active_signaling_network,lr_network = lr_network, sig_network = sig_network, gr_network = gr_network)
head(data_source_network) 
## # A tibble: 6 × 5
##   from  to     source               database       layer     
##   <chr> <chr>  <chr>                <chr>          <chr>     
## 1 ABCA1 CRK    harmonizome_GEO_GENE harmonizome_gr regulatory
## 2 ABCA1 DGKZ   harmonizome_GEO_GENE harmonizome_gr regulatory
## 3 ABCA1 INPP5D harmonizome_GEO_GENE harmonizome_gr regulatory
## 4 ACHE  CRK    harmonizome_GEO_GENE harmonizome_gr regulatory
## 5 ACHE  PIK3CA harmonizome_GEO_GENE harmonizome_gr regulatory
## 6 ACHE  PTPN11 harmonizome_GEO_GENE harmonizome_gr regulatory

export and save network and metadata:

output_path = "Results/all_NK_cells/"
dir.create(output_path, recursive = TRUE)
## Warning in dir.create(output_path, recursive = TRUE): 'Results/all_NK_cells'
## already exists
write_output = TRUE # change to TRUE for writing output
# weighted networks ('import network' in Cytoscape)
save_nets <- function(active_signaling_network, ligands, targets, output_path, write_output){
  if(write_output){
    bind_rows(active_signaling_network$sig %>% mutate(layer = "signaling"), active_signaling_network$gr %>% mutate(layer = "regulatory")) %>% write_tsv(paste0(output_path,"weighted_signaling_network.txt")) 
  }
  
  # networks with information of supporting data sources ('import network' in Cytoscape)
  if(write_output){
  data_source_network %>% write_tsv(paste0(output_path,"data_source_network.txt"))
  }
  
  # Node annotation table ('import table' in Cytoscape)
  specific_annotation_tbl = bind_rows(
    tibble(gene = ligands, annotation = "ligand"),
    tibble(gene = targets, annotation = "target"),
    tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(c(targets,ligands)) %>% intersect(lr_network$to %>% unique()), annotation = "receptor"),
    tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(c(targets,ligands)) %>% intersect(gr_network$from %>% unique()) %>% setdiff(c(data_source_network$from, data_source_network$to) %>% unique() %>% intersect(lr_network$to %>% unique())),annotation = "transcriptional regulator")
  )
  non_specific_annotation_tbl = tibble(gene = c(data_source_network$from, data_source_network$to) %>% unique() %>% setdiff(specific_annotation_tbl$gene), annotation = "signaling mediator")
  
  if(write_output){
  bind_rows(specific_annotation_tbl,non_specific_annotation_tbl) %>% write_tsv(paste0(output_path,"annotation_table.txt"))
  }
}

Sub-networks for each ligand:

sub_nets <- function(ligand_tf_matrix,ligand, targets){
  
  
  active_signaling_network = get_ligand_signaling_path(ligand_tf_matrix = ligand_tf_matrix,
                                                       ligands_all = ligand,
                                                       targets_all = targets,
                                                       weighted_networks = weighted_networks)
  
  data_source_network = infer_supporting_datasources(signaling_graph_list = active_signaling_network,
                                                     lr_network = lr_network,
                                                     sig_network = sig_network,
                                                     gr_network = gr_network)



  active_signaling_network_min_max = active_signaling_network
  active_signaling_network_min_max$sig = active_signaling_network_min_max$sig %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
  active_signaling_network_min_max$gr = active_signaling_network_min_max$gr %>% mutate(weight = ((weight-min(weight))/(max(weight)-min(weight))) + 0.75)
  
  graph_min_max = diagrammer_format_signaling_graph(signaling_graph_list = active_signaling_network_min_max, ligands_all = ligand, targets_all = targets, sig_color = "indianred", gr_color = "steelblue")
  
  # To render the graph: uncomment following line of code
  graph = DiagrammeR::render_graph(graph_min_max, layout = "kk", title = ligand)

  output_path = paste0("Results/all_NK_cells/", ligand, "/")
  dir.create(output_path, recursive = TRUE)
  write_output = TRUE # change to TRUE for writing output
  
  save_nets(active_signaling_network,ligands = ligand, targets = targets, output_path, write_output)
  return(graph)
}
plot_list <- htmltools::tagList()

for (l in 1:length(ligands)){
  graph = sub_nets(ligand_tf_matrix, ligands[l], geneset_oi)
  plot_list[[l]] <- graph
}
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the DiagrammeR package.
##   Please report the issue at
##   <]8;;https://github.com/rich-iannone/DiagrammeR/issueshttps://github.com/rich-iannone/DiagrammeR/issues]8;;>.
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/B2M' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CALM1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CCN2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD22' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD40LG' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD48' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD55' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CD99' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CDH1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CEACAM1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CLEC2B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CLEC2D' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/CXCL12' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/FAM3C' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/GNAS' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/GZMB' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-A' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-C' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DPB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQA1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQA2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DQB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRA' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-DRB5' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HLA-E' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/HMGB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/IL1B' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS3' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LGALS9' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LILRB4' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/LTF' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/MIF' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NECTIN1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NECTIN2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/NRG1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/SIRPA' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TGFB1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TGM2' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/THBS1' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TNFSF13' already exists
## Warning in dir.create(output_path, recursive = TRUE):
## 'Results/all_NK_cells/TNFSF13B' already exists
plot_list

Conclusions

The interactions with LAG3, TIGIT, LAIR1 and other ICR looks interesting. Because of problems visualizing the whole network on R I’m going to continue the analysis on python.

I built the sub-networks with all the target genes (all the genes included in those regulons that were including ICRs mediators). This was more for testing and visualization purposes. Selecting specific ligands or, idndirectly receptors, and some target genes (e.g. SHP1/SHP2), maybe we could identify smaller sub-networks with clearer pathways. Now im going to focus on the network dedicated to only exh NK cells in MM and check which ligands and target genes could better suit our aims.